########################
### Replication File for "AFFILIATION BIAS IN ARBITRATION: AN EXPERIMENTAL APPROACH"
### Puig, Strezhnev
########################
### Note: Be sure to set your working directory to the location of the three .csv data files

### Clear environment
rm(list=ls())

##### Libraries
library(Hmisc)
library(gtools)
library(gplots)
library(ggplot2)

#### Globals
nIter <- 5000
random.seed <- 1628 # Generated via "Random.org" [1-5000]

############################################
##### First Experiment
############################################

#################################
## Data pre-processing - Experiment 1
#################################

## Load in cleaned results file
data <- read.csv("experiment1_public.csv")

##### Note: -99 = missing

## How many initial respondents
sum(!is.na(data$ConsentForm)&data$ConsentForm!=2)

### Drop anyone who didn't answer 1 (Yes) to the consent form
data <- subset(data, !is.na(data$ConsentForm)&data$ConsentForm!=2)

## Make respondent ID a character
data$ResponseID <- as.character(data$ResponseID)

## Merge reimbursement questions
## Investor-State
data$isClaimantReimbursesAll[!is.na(data$isClaimantReimbursesAll2)] <- data$isClaimantReimbursesAll2[!is.na(data$isClaimantReimbursesAll2)]
data$isClaimantReimbursesSome[!is.na(data$isClaimantReimbursesSome2)] <- data$isClaimantReimbursesSome2[!is.na(data$isClaimantReimbursesSome2)]
data$isSplitExpenses[!is.na(data$isSplitExpenses2)] <- data$isSplitExpenses2[!is.na(data$isSplitExpenses2)]
data$isRespondentReimbursesSome[!is.na(data$isRespondentReimbursesSome2)] <- data$isRespondentReimbursesSome2[!is.na(data$isRespondentReimbursesSome2)]
data$isRespondentReimbursesAll[!is.na(data$isRespondentReimbursesAll2)] <- data$isRespondentReimbursesAll2[!is.na(data$isRespondentReimbursesAll2)]

## Recode outcome levels to be more legible
levels(data$isOutcome)[grepl("respondent did NOT unfairly treat", levels(data$isOutcome))] <- "Respondent wins"
levels(data$isOutcome)[grepl("claims were manifestly without legal merit", levels(data$isOutcome))] <- "Claimant without legal merit"
levels(data$isOutcome)[grepl("present dispute was outside of the jurisdiction", levels(data$isOutcome))] <- "Outside jurisdiction"
levels(data$isOutcome)[grepl("respondent unfairly treated and wrongfully expropriated", levels(data$isOutcome))] <- "Claimant wins"

## Recode blank level as "blind appointment"
levels(data$isAppointer)[grepl("claimant", levels(data$isAppointer))] <- "Claimant"
levels(data$isAppointer)[grepl("respondent", levels(data$isAppointer))] <- "Respondent"
levels(data$isAppointer)[grepl("parties", levels(data$isAppointer))] <- "Parties"
levels(data$isAppointer)[grepl(".", levels(data$isAppointer), fixed=T)] <- "Blind"

####### Recode -99 as NA
data$genderMale[data$genderMale == -99] <- NA
data$genderFemale[data$genderFemale == -99] <- NA

####### Recode cost share variable 
data$isReimbursement <- NA
data$isReimbursement[data$isClaimantReimbursesAll == 1] <- 1
data$isReimbursement[data$isClaimantReimbursesSome == 1] <- 2
data$isReimbursement[data$isSplitExpenses == 1] <- 3
data$isReimbursement[data$isRespondentReimbursesSome == 1] <- 4
data$isReimbursement[data$isRespondentReimbursesAll == 1] <- 5

#####################################################
### Covariates of respondents

### Gender
data$genderMale[data$genderMale == -99] <- NA

### Employment
data$employmentPrivate[data$employmentPrivate == -99] <- NA
data$employmentProfessor[data$employmentProfessor == -99] <- NA
data$employmentGovt[data$employmentGovt == -99] <- NA
data$employmentOther[data$employmentOther== -99] <- NA
data$employmentCoded[data$employmentCoded == ""] <- NA

### Served on IS dispute?
data$isDisputeYes[data$isDisputeYes == -99] <- NA

######################################
### Treatment variables
######################################

### Who won the dispute?
data$isWinner <- as.factor(ifelse(data$isOutcome == "Claimant wins", "Claimant", "Respondent"))
data$isWinner <- relevel(data$isWinner, "Respondent")

### Were the parties appointed by the winner vs. loser?
data$isAppointerSide <- as.character(data$isAppointer)
data$isAppointerSide[(data$isAppointer == "Claimant"& data$isWinner == "Claimant")|
                         (data$isAppointer == "Respondent"& data$isWinner == "Respondent")] <- "Appointed by Winner"
data$isAppointerSide[(data$isAppointer == "Claimant"& data$isWinner == "Respondent")|
                         (data$isAppointer == "Respondent"& data$isWinner == "Claimant")] <- "Appointed by Loser"
data$isAppointerSide <- as.factor(data$isAppointerSide)

data$isAppointerSide <- relevel(data$isAppointerSide, "Appointed by Loser")

### Reimbursement - reframe in terms of the winner of the dispute
data$isReimbursementFolded <- data$isReimbursement
data$isReimbursementFolded[data$isWinner == "Respondent"] <- 6 - data$isReimbursement[data$isWinner == "Respondent"]

###### Subset to units that finished the survey
completed.data <- subset(data, Finished == 1)

############################################
#### Effect of appointer
############################################

### Subset to cases where respondents answered the outcome question on costs
completed.data.full <- subset(completed.data, !is.na(completed.data$isReimbursementFolded))

### Number of observations in total
nrow(completed.data.full)

### Check for whether randomization holds, p > .05 in all cases.
chisq.test(table(completed.data.full$isAppointer), p=c(1/4, 1/4, 1/4, 1/4)) ## Appointer
chisq.test(table(completed.data.full$isOutcome), p=c(1/4, 1/4, 1/4, 1/4)) ## Outcome
chisq.test(table(completed.data.full$isClaimant), p=c(1/2, 1/2)) ## Claimant
chisq.test(table(completed.data.full$isRespondent), p=c(1/2, 1/2)) ## Respondent

### Parse outcome into a bunch of binary indicators
completed.data.full$isSplit <- as.integer(completed.data.full$isReimbursementFolded == 3)
completed.data.full$isSome <- as.integer(completed.data.full$isReimbursementFolded == 4)
completed.data.full$isAll <- as.integer(completed.data.full$isReimbursementFolded == 5)
completed.data.full$isAny <- as.integer(completed.data.full$isReimbursementFolded == 4 | completed.data.full$isReimbursementFolded == 5)

######################
#### Estimate treament effects and bootstrap the standard errors
######################

#### Placeholder for holding point estimates and standard errors

results_frame_exp1 <- data.frame(treatment=c(), outcome=c(), estimate=c(), standard_error = c(), pval = c())

######### Effect of appointment by winner vs. appointment by loser

### Set the random seed before sampling
set.seed(random.seed)

## Matrix to store marginal effect
appointerWinner <- matrix(nrow=nIter, ncol=4)

## Make "Appointed by Loser" the baseline level
completed.data.full$isAppointerSide <- relevel(completed.data.full$isAppointerSide, "Appointed by Loser")

# Bootstrap routine
for (iter in 1:nIter){
  
  # Resample rows of data
  if (iter != 1){ 
    indices <- sample(1:nrow(completed.data.full), nrow(completed.data.full), replace=T) 
    boot.data <- completed.data.full[indices,]
  }else{ # Keep the data the same for the first bootstrap iteration
    boot.data <- completed.data.full
  }
  
  # Difference-in-means (fully-saturated regression models)
  boot_split <- lm(isSplit ~ isAppointerSide, data=boot.data) ## Split costs
  boot_some <- lm(isSome ~  isAppointerSide, data=boot.data) ## Some costs
  boot_all <- lm(isAll ~  isAppointerSide, data=boot.data) ## All Costs
  boot_any <- lm(isAny ~ isAppointerSide, data=boot.data) ## All or Some Costs
  
  ## Store coefficients denoting differences-in-means
  appointerWinner[iter,] <- c(coefficients(boot_split)["isAppointerSideAppointed by Winner"], coefficients(boot_some)["isAppointerSideAppointed by Winner"], coefficients(boot_all)["isAppointerSideAppointed by Winner"],coefficients(boot_any)["isAppointerSideAppointed by Winner"])

}

## Point Estimate
appointerWinner_point <- appointerWinner[1,]
## SD of bootstraps is our estimated SE
appointerWinner_se <- apply(appointerWinner, 2, function(x) sd(x))
## CIs are quantiles of the bootstrap distribution
appointerWinner_ci <- apply(appointerWinner, 2, function(x) quantile(x, c(.025, .975)))
appointerWinner_ci90 <- apply(appointerWinner, 2, function(x) quantile(x, c(.05, .95)))
## P-value from SE estimate of bootstrap (sampling distribution is asymptotically normal, so this doesn't differ much from just taking quantiles)
appointerWinner_pval <- sapply(appointerWinner_point/appointerWinner_se, function(x) 2*pnorm(-abs(x)))

### Output point/SE/CI estimates
round(appointerWinner_point, 3)
round(appointerWinner_se, 3)
round(appointerWinner_pval, 3)

### Store in regression tables

results_frame_exp1 <- rbind(results_frame_exp1,
                       data.frame(treatment=c("Appointed by Winner v. Loser","Appointed by Winner v. Loser","Appointed by Winner v. Loser"), outcome=c("Split Costs", "Some Costs", "All Costs"), estimate=round(appointerWinner_point, 3)[1:3], 
                                  standard_error = round(appointerWinner_se, 3)[1:3], pval = round(appointerWinner_pval, 3)[1:3]))

### Sample sizes for each treatment condition
table(completed.data.full$isAppointerSide)

########################################
#### Plot of marginal effects for Winner vs. Loser case
########################################

### Put estimates in a data frame
winnerloser_effect <- data.frame(choice = c( "Parties Pay Own Costs","Loser Reimburses Some","Loser Reimburses All"),
                                 point= appointerWinner[1,c(1:3)],
                                 lower95 = apply(appointerWinner, 2, function(x) quantile(x, c(.025, .975)))[1,c(1:3)], upper95=
                                   apply(appointerWinner, 2, function(x) quantile(x, c(.025, .975)))[2,c(1:3)],
                                 lower90 = apply(appointerWinner, 2, function(x) quantile(x, c(.05, .95)))[1,c(1:3)], upper90=
                                   apply(appointerWinner, 2, function(x) quantile(x, c(.05, .95)))[2,c(1:3)])

### Plot output
png("appointed_by_winner_exp1.png", width=900, height=350)
p = ggplot(winnerloser_effect, aes(y=point, x=choice))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower90,ymax=upper90), size= 4)
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Change in the probability respondents select the outcome", limits=c(-.4, .4))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + ggtitle("Appointed by winning party versus losing party")
p = p + theme_bw1()
print(p)
dev.off()
ggsave("appointed_by_winner_exp1.pdf",width=12, height=4)

#################
## Blind v. non-blind
###################

######### Bootstrap

## Set random seed
set.seed(random.seed)

## Matrix to store marginal effect
appointerParties <- matrix(nrow=nIter, ncol=4)

## New baseline level "Blind"
completed.data.full$isAppointerSide <- relevel(completed.data.full$isAppointerSide, "Blind")

# Bootstrap routine
for (iter in 1:nIter){

  # Resample rows of data
  if (iter != 1){ # Keep the data the same for the first estimate
    indices <- sample(1:nrow(completed.data.full), nrow(completed.data.full), replace=T) 
    boot.data <- completed.data.full[indices,]
  }else{
    boot.data <- completed.data.full
  }
  
  # Difference-in-means
  boot_split <- lm(isSplit ~ isAppointerSide, data=boot.data)
  boot_some <- lm(isSome ~  isAppointerSide, data=boot.data)
  boot_all <- lm(isAll ~  isAppointerSide, data=boot.data)
  boot_any <- lm(isAny ~ isAppointerSide, data=boot.data)
  
  appointerParties[iter,] <- c(coefficients(boot_split)["isAppointerSideParties"], coefficients(boot_some)["isAppointerSideParties"], coefficients(boot_all)["isAppointerSideParties"],coefficients(boot_any)["isAppointerSideParties"])
  
}

appointerParties_point <- appointerParties[1,]
appointerParties_se <- apply(appointerParties, 2, function(x) sd(x))
appointerParties_ci <- apply(appointerParties, 2, function(x) quantile(x, c(.025, .975)))
appointerParties_ci90 <- apply(appointerParties, 2, function(x) quantile(x, c(.05, .95)))
appointerParties_pval <- sapply(appointerParties_point/appointerParties_se, function(x) 2*pnorm(-abs(x)))

round(appointerParties_point, 3)
round(appointerParties_se, 3)
round(appointerParties_pval, 3)

### Store points
results_frame_exp1 <- rbind(results_frame_exp1,
                            data.frame(treatment=c("Appointed by Parties v. Blind","Appointed by Parties v. Blind","Appointed by Parties v. Blind"), outcome=c("Split Costs", "Some Costs", "All Costs"), estimate=round(appointerParties_point, 3)[1:3], 
                                       standard_error = round(appointerParties_se, 3)[1:3], pval = round(appointerParties_pval, 3)[1:3]))


### Sample sizes for each treatment condition
table(completed.data.full$isAppointerSide)

########################################
#### Plot of marginal effects for Winner vs. Loser case
########################################


partyblind_effect <- data.frame(choice = c( "Parties Pay Own Costs","Loser Reimburses Some","Loser Reimburses All"),
                                 point= appointerParties[1,c(1:3)],
                                 lower95 = apply(appointerParties, 2, function(x) quantile(x, c(.025, .975)))[1,c(1:3)], upper95=
                                   apply(appointerParties, 2, function(x) quantile(x, c(.025, .975)))[2,c(1:3)],
                                 lower90 = apply(appointerParties, 2, function(x) quantile(x, c(.05, .95)))[1,c(1:3)], upper90=
                                   apply(appointerParties, 2, function(x) quantile(x, c(.05, .95)))[2,c(1:3)])



png("appointed_blind_exp1.png", width=900, height=350)
p = ggplot(partyblind_effect, aes(y=point, x=choice))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower90,ymax=upper90), size= 4)
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Change in the probability respondents select the outcome", limits=c(-.4, .4))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + ggtitle("Appointed by the parties versus blind")
p = p + theme_bw1()
print(p)
dev.off()
ggsave("appointed_blind_exp1.pdf",width=12, height=4)

########################################## Print results in a table
## Results table - Experiment 1
results_frame_exp1

#######################################################################################
###### Compare sampled population to the target Population
#######################################################################################

####### Complete cases (rename the data)
folded.data <- completed.data.full

### Load dataset on the targeted arbitrators
target_arbitrators <- read.csv("icsid_2010-2016_arbitrators.csv")

### Rename variables/make comparable with our primary dataset
target_arbitrators$employmentCoded <- target_arbitrators$Profession.Primary
folded.data$Gender <- as.factor(ifelse(folded.data$genderMale, "M", "F"))
folded.data$employmentCoded <- factor(folded.data$employmentCoded)

target_arbitrators$legalOrigin <- NA
target_arbitrators$legalOrigin[target_arbitrators$anyCommon == 1 & target_arbitrators$anyCivil == 0] <- "Common Law"
target_arbitrators$legalOrigin[target_arbitrators$anyCommon == 0 & target_arbitrators$anyCivil == 1] <- "Civil Law"
target_arbitrators$legalOrigin[target_arbitrators$anyCommon == 1 & target_arbitrators$anyCivil == 1] <- "Both"

folded.data$legalOrigin <- NA
folded.data$legalOrigin[folded.data$anyCommon == 1 & folded.data$anyCivil == 0] <- "Common Law"
folded.data$legalOrigin[folded.data$anyCommon == 0 & folded.data$anyCivil == 1] <- "Civil Law"
folded.data$legalOrigin[folded.data$anyCommon == 1 & folded.data$anyCivil == 1] <- "Both"

folded.data$legalOrigin <- as.factor(folded.data$legalOrigin)

### Subset data to just the observations that had responses to all of the demographic questions.
folded.data.complete <- subset(folded.data, !is.na(folded.data$legalOrigin)&!is.na(folded.data$Gender)&!is.na(folded.data$employmentCoded))

### What are the population proportions in the three variables we care about
targetList <- list()
targetList[["Gender"]] <- table(target_arbitrators$Gender)/sum(table(target_arbitrators$Gender))
targetList[["legalOrigin"]] <- table(target_arbitrators$legalOrigin)/sum(table(target_arbitrators$legalOrigin))
targetList[["employmentCoded"]] <- table(target_arbitrators$employmentCoded)/sum(table(target_arbitrators$employmentCoded))

#####
### Gender 
####

#################### Distribution of gender

gender_freq <- data.frame(category = c("Female", "Male"),
                          point = as.integer(table(folded.data.complete$Gender)))

gender_freq <- cbind(gender_freq, binconf(x = gender_freq$point, n=sum(gender_freq$point)))

png("gender_frequencies_newsamp.png", width=400, height=350)
par(cex=1)
par(cex.main=1.3)
par(las=2)
par(mar=c(5,6,4,2)+0.4)
bplot <- barplot2(gender_freq$PointEst, horiz=T, width=.2,
                  names.arg=c(paste("Female\n(", gender_freq$point[1], " total)", sep=""), 
                              paste("Male\n(", gender_freq$point[2], " total)", sep="")), xlab = "Share of responses", lwd=2,
                  main=paste("Gender of respondents in sample, n=", sum(gender_freq$point), sep=""), xlim=c(0,1), plot.ci=T, ci.lwd=2, ci.width=.8, ci.color="grey20",
                  ci.l = gender_freq$Lower, 
                  ci.u = gender_freq$Upper )
points(y=bplot+.1, x=as.vector(table(target_arbitrators$Gender)/sum(table(target_arbitrators$Gender))), cex=2, pch=4, lwd=4)
segments(x0 = c(0,0), y0=bplot+.1, x1= as.vector(table(target_arbitrators$Gender)/sum(table(target_arbitrators$Gender))), y1 = bplot+.1, lty=2, lwd=4)
dev.off()
### PDF version
pdf("gender_frequencies_newsamp.pdf", width=5.4,height=4)
par(cex=1)
par(cex.main=1.3)
par(las=2)
par(mar=c(5,6,4,2)+0.4)
bplot <- barplot2(gender_freq$PointEst, horiz=T, width=.2,
                  names.arg=c(paste("Female\n(", gender_freq$point[1], " total)", sep=""), 
                              paste("Male\n(", gender_freq$point[2], " total)", sep="")), xlab = "Share of responses", lwd=2,
                  main=paste("Gender of respondents in sample, n=", sum(gender_freq$point), sep=""), xlim=c(0,1), plot.ci=T, ci.lwd=2, ci.width=.8, ci.color="grey20",
                  ci.l = gender_freq$Lower, 
                  ci.u = gender_freq$Upper )
points(y=bplot+.1, x=as.vector(table(target_arbitrators$Gender)/sum(table(target_arbitrators$Gender))), cex=2, pch=4, lwd=4)
segments(x0 = c(0,0), y0=bplot+.1, x1= as.vector(table(target_arbitrators$Gender)/sum(table(target_arbitrators$Gender))), y1 = bplot+.1, lty=2, lwd=4)
dev.off()

#############
### Is there a gender interaction? Winner v. Loser
###############

## Set random seed
set.seed(random.seed)

### Dataset with complete responses
folded.data.complete.resp <- subset(folded.data.complete, !is.na(completed.data$isReimbursementFolded))

## Matrix to store interaction effect
appointerWinner_Gender <- matrix(nrow=nIter, ncol=4)

## Make "Appointed by Loser" the baseline level
folded.data.complete.resp$isAppointerSide <- relevel(folded.data.complete.resp$isAppointerSide, "Appointed by Loser")

# Bootstrap routine
for (iter in 1:nIter){
  
  # Resample rows of data
  if (iter != 1){ 
    indices <- sample(1:nrow(folded.data.complete.resp), nrow(folded.data.complete.resp), replace=T) 
    boot.data <- folded.data.complete.resp[indices,]
  }else{ # Keep the data the same for the first bootstrap iteration
    boot.data <- folded.data.complete.resp
  }
  
  # Difference-in-means (fully-saturated regression models)
  boot_split <- lm(isSplit ~ isAppointerSide*Gender, data=boot.data) ## Split costs
  boot_some <- lm(isSome ~  isAppointerSide*Gender, data=boot.data) ## Some costs
  boot_all <- lm(isAll ~  isAppointerSide*Gender, data=boot.data) ## All Costs
  boot_any <- lm(isAny ~ isAppointerSide*Gender, data=boot.data) ## All or Some Costs
  
  ## Store coefficients denoting difference-in-difference in means between men and women
  appointerWinner_Gender[iter,] <- c(coefficients(boot_split)["isAppointerSideAppointed by Winner:GenderM"], coefficients(boot_some)["isAppointerSideAppointed by Winner:GenderM"], coefficients(boot_all)["isAppointerSideAppointed by Winner:GenderM"],coefficients(boot_any)["isAppointerSideAppointed by Winner:GenderM"])
  
}

## Point Estimate
appointerWinner_Gender_point <- appointerWinner_Gender[1,]
## SD of bootstraps is our estimated SE
appointerWinner_Gender_se <- apply(appointerWinner_Gender, 2, function(x) sd(x, na.rm=T))
## CIs are quantiles of the bootstrap distribution
appointerWinner_Gender_ci <- apply(appointerWinner_Gender, 2, function(x) quantile(x, c(.025, .975), na.rm=T))
appointerWinner_Gender_ci90 <- apply(appointerWinner_Gender, 2, function(x) quantile(x, c(.05, .95), na.rm=T))
## P-value from SE estimate of bootstrap (sampling distribution is asymptotically normal, so this doesn't differ much from just taking quantiles)
appointerWinner_Gender_pval <- sapply(appointerWinner_Gender_point/appointerWinner_Gender_se, function(x) 2*pnorm(-abs(x)))

## No Statistically significant interaction with gender
appointerWinner_Gender_pval

################### Frequencies of sample Legal origins

legal_freq <- data.frame(category = c("Both", "Civil Law", "Common Law"),
                         point = as.integer(table(folded.data.complete$legalOrigin)))
legal_freq <- cbind(legal_freq, binconf(x = legal_freq$point, n=sum(legal_freq$point)))

png("legal_frequencies_newsamp.png", width=400, height=350)
par(cex=1)
par(cex.main=1.3)
par(las=2)
par(mar=c(5,6,4,2)+0.4)
bplot <- barplot2(legal_freq$PointEst, horiz=T, width=.2,
                  names.arg=c(paste("Both\n(", legal_freq$point[1], " total)", sep=""), 
                              paste("Civil\nLaw\n(", legal_freq$point[2], " total)", sep=""),
                              paste("Common\nLaw\n(", legal_freq$point[3], " total)", sep="")), xlab = "Share of responses", lwd=2,
                  main=paste("Legal origins of respondents\nin sample, n=", sum(legal_freq$point), sep=""), xlim=c(0,1), plot.ci=T, ci.lwd=2, ci.width=.8, ci.color="grey20",
                  ci.l = legal_freq$Lower, 
                  ci.u = legal_freq$Upper )
points(y=bplot+.1, x=as.vector(table(target_arbitrators$legalOrigin)/sum(table(target_arbitrators$legalOrigin))), cex=2, pch=4, lwd=4)
segments(x0 = c(0,0), y0=bplot+.1, x1= as.vector(table(target_arbitrators$legalOrigin)/sum(table(target_arbitrators$legalOrigin))), y1 = bplot+.1, lty=2, lwd=4)
dev.off()
### PDF Version
pdf("legal_frequencies_newsamp.pdf", width=5.4, height=4)
par(cex=1)
par(cex.main=1.3)
par(las=2)
par(mar=c(5,6,4,2)+0.4)
bplot <- barplot2(legal_freq$PointEst, horiz=T, width=.2,
                  names.arg=c(paste("Both\n(", legal_freq$point[1], " total)", sep=""), 
                              paste("Civil\nLaw\n(", legal_freq$point[2], " total)", sep=""),
                              paste("Common\nLaw\n(", legal_freq$point[3], " total)", sep="")), xlab = "Share of responses", lwd=2,
                  main=paste("Legal origins of respondents\nin sample, n=", sum(legal_freq$point), sep=""), xlim=c(0,1), plot.ci=T, ci.lwd=2, ci.width=.8, ci.color="grey20",
                  ci.l = legal_freq$Lower, 
                  ci.u = legal_freq$Upper )
points(y=bplot+.1, x=as.vector(table(target_arbitrators$legalOrigin)/sum(table(target_arbitrators$legalOrigin))), cex=2, pch=4, lwd=4)
segments(x0 = c(0,0), y0=bplot+.1, x1= as.vector(table(target_arbitrators$legalOrigin)/sum(table(target_arbitrators$legalOrigin))), y1 = bplot+.1, lty=2, lwd=4)
dev.off()

###############
### Is there a Common Law interaction? Winner v. Loser
###############

## Set random seed
set.seed(random.seed)

### Dataset with complete responses
folded.data.complete.resp <- subset(folded.data.complete, !is.na(completed.data$isReimbursementFolded))

## Matrix to store interaction effect
appointerWinner_CommonLaw <- matrix(nrow=nIter, ncol=4)

## Make "Appointed by Loser" the baseline level
folded.data.complete.resp$isAppointerSide <- relevel(folded.data.complete.resp$isAppointerSide, "Appointed by Loser")

# Bootstrap routine
for (iter in 1:nIter){
  
  # Resample rows of data
  if (iter != 1){ 
    indices <- sample(1:nrow(folded.data.complete.resp), nrow(folded.data.complete.resp), replace=T) 
    boot.data <- folded.data.complete.resp[indices,]
  }else{ # Keep the data the same for the first bootstrap iteration
    boot.data <- folded.data.complete.resp
  }
  
  # Difference-in-means (fully-saturated regression models)
  boot_split <- lm(isSplit ~ isAppointerSide*I(legalOrigin == "Common Law"), data=boot.data) ## Split costs
  boot_some <- lm(isSome ~  isAppointerSide*I(legalOrigin == "Common Law"), data=boot.data) ## Some costs
  boot_all <- lm(isAll ~  isAppointerSide*I(legalOrigin == "Common Law"), data=boot.data) ## All Costs
  boot_any <- lm(isAny ~ isAppointerSide*I(legalOrigin == "Common Law"), data=boot.data) ## All or Some Costs
  
  ## Store coefficients denoting difference-in-difference in means between men and women
  appointerWinner_CommonLaw[iter,] <- c(coefficients(boot_split)['isAppointerSideAppointed by Winner:I(legalOrigin == "Common Law")TRUE'], 
                                        coefficients(boot_some)['isAppointerSideAppointed by Winner:I(legalOrigin == "Common Law")TRUE'], 
                                        coefficients(boot_all)['isAppointerSideAppointed by Winner:I(legalOrigin == "Common Law")TRUE'],
                                        coefficients(boot_any)['isAppointerSideAppointed by Winner:I(legalOrigin == "Common Law")TRUE'])
  
}

## Point Estimate
appointerWinner_CommonLaw_point <- appointerWinner_CommonLaw[1,]
## SD of bootstraps is our estimated SE
appointerWinner_CommonLaw_se <- apply(appointerWinner_CommonLaw, 2, function(x) sd(x, na.rm=T))
## CIs are quantiles of the bootstrap distribution
appointerWinner_CommonLaw_ci <- apply(appointerWinner_CommonLaw, 2, function(x) quantile(x, c(.025, .975), na.rm=T))
appointerWinner_CommonLaw_ci90 <- apply(appointerWinner_CommonLaw, 2, function(x) quantile(x, c(.05, .95), na.rm=T))
## P-value from SE estimate of bootstrap (sampling distribution is asymptotically normal, so this doesn't differ much from just taking quantiles)
appointerWinner_CommonLaw_pval <- sapply(appointerWinner_CommonLaw_point/appointerWinner_CommonLaw_se, function(x) 2*pnorm(-abs(x)))

## No Statistically significant interaction with Legal origins
appointerWinner_CommonLaw_pval


################### Frequencies of sample professions/backgrounds

background_freq <- data.frame(category = c("Academic", "Private", "Public"),
                              point = as.integer(table(folded.data.complete$employmentCoded)))
background_freq <- cbind(background_freq, binconf(x = background_freq$point, n=sum(background_freq$point)))

png("background_frequencies_newsamp.png", width=400, height=350)
par(cex=1)
par(cex.main=1.3)
par(las=2)
par(mar=c(5,6,4,2)+0.4)
bplot <- barplot2(background_freq$PointEst, horiz=T, width=.2,
                  names.arg=c(paste("Academic\n(", background_freq$point[1], " total)", sep=""), 
                              paste("Private\n(", background_freq$point[2], " total)", sep=""),
                              paste("Public\n(", background_freq$point[3], " total)", sep="")), xlab = "Share of responses", lwd=2,
                  main=paste("Employment of respondents\nin sample, n=", sum(background_freq$point), sep=""), xlim=c(0,1), plot.ci=T, ci.lwd=2, ci.width=.8, ci.color="grey20",
                  ci.l = background_freq$Lower, 
                  ci.u = background_freq$Upper )
points(y=bplot+.1, x=as.vector(table(target_arbitrators$employmentCoded)/sum(table(target_arbitrators$employmentCoded))), cex=2, pch=4, lwd=4)
segments(x0 = c(0,0), y0=bplot+.1, x1= as.vector(table(target_arbitrators$employmentCoded)/sum(table(target_arbitrators$employmentCoded))), y1 = bplot+.1, lty=2, lwd=4)
dev.off()
### PDF version
pdf("background_frequencies_newsamp.pdf", width=5.4, height=4)
par(cex=1)
par(cex.main=1.3)
par(las=2)
par(mar=c(5,6,4,2)+0.4)
bplot <- barplot2(background_freq$PointEst, horiz=T, width=.2,
                  names.arg=c(paste("Academic\n(", background_freq$point[1], " total)", sep=""), 
                              paste("Private\n(", background_freq$point[2], " total)", sep=""),
                              paste("Public\n(", background_freq$point[3], " total)", sep="")), xlab = "Share of responses", lwd=2,
                  main=paste("Employment of respondents\nin sample, n=", sum(background_freq$point), sep=""), xlim=c(0,1), plot.ci=T, ci.lwd=2, ci.width=.8, ci.color="grey20",
                  ci.l = background_freq$Lower, 
                  ci.u = background_freq$Upper )
points(y=bplot+.1, x=as.vector(table(target_arbitrators$employmentCoded)/sum(table(target_arbitrators$employmentCoded))), cex=2, pch=4, lwd=4)
segments(x0 = c(0,0), y0=bplot+.1, x1= as.vector(table(target_arbitrators$employmentCoded)/sum(table(target_arbitrators$employmentCoded))), y1 = bplot+.1, lty=2, lwd=4)
dev.off()

###############
### Is there a Private Sector interaction? Winner v. Loser
###############

## Set random seed
set.seed(random.seed)

### Dataset with complete responses
folded.data.complete.resp <- subset(folded.data.complete, !is.na(completed.data$isReimbursementFolded))

## Matrix to store interaction effect
appointerWinner_Private <- matrix(nrow=nIter, ncol=4)

## Make "Appointed by Loser" the baseline level
folded.data.complete.resp$isAppointerSide <- relevel(folded.data.complete.resp$isAppointerSide, "Appointed by Loser")

# Bootstrap routine
for (iter in 1:nIter){
  
  # Resample rows of data
  if (iter != 1){ 
    indices <- sample(1:nrow(folded.data.complete.resp), nrow(folded.data.complete.resp), replace=T) 
    boot.data <- folded.data.complete.resp[indices,]
  }else{ # Keep the data the same for the first bootstrap iteration
    boot.data <- folded.data.complete.resp
  }
  
  # Difference-in-means (fully-saturated regression models)
  boot_split <- lm(isSplit ~ isAppointerSide*I(employmentCoded == "Private"), data=boot.data) ## Split costs
  boot_some <- lm(isSome ~  isAppointerSide*I(employmentCoded == "Private"), data=boot.data) ## Some costs
  boot_all <- lm(isAll ~  isAppointerSide*I(employmentCoded == "Private"), data=boot.data) ## All Costs
  boot_any <- lm(isAny ~ isAppointerSide*I(employmentCoded == "Private"), data=boot.data) ## All or Some Costs
  
  ## Store coefficients denoting difference-in-difference in means between men and women
  appointerWinner_Private[iter,] <- c(coefficients(boot_split)['isAppointerSideAppointed by Winner:I(employmentCoded == "Private")TRUE'], 
                                        coefficients(boot_some)['isAppointerSideAppointed by Winner:I(employmentCoded == "Private")TRUE'], 
                                        coefficients(boot_all)['isAppointerSideAppointed by Winner:I(employmentCoded == "Private")TRUE'],
                                        coefficients(boot_any)['isAppointerSideAppointed by Winner:I(employmentCoded == "Private")TRUE'])
  
}

## Point Estimate
appointerWinner_Private_point <- appointerWinner_Private[1,]
## SD of bootstraps is our estimated SE
appointerWinner_Private_se <- apply(appointerWinner_Private, 2, function(x) sd(x, na.rm=T))
## CIs are quantiles of the bootstrap distribution
appointerWinner_Private_ci <- apply(appointerWinner_Private, 2, function(x) quantile(x, c(.025, .975), na.rm=T))
appointerWinner_Private_ci90 <- apply(appointerWinner_Private, 2, function(x) quantile(x, c(.05, .95), na.rm=T))
## P-value from SE estimate of bootstrap (sampling distribution is asymptotically normal, so this doesn't differ much from just taking quantiles)
appointerWinner_Private_pval <- sapply(appointerWinner_Private_point/appointerWinner_Private_se, function(x) 2*pnorm(-abs(x)))

## No Statistically significant interaction with Private sector employment
appointerWinner_Private_pval

### END EXPERIMENT 1

################################################################
###############################################################
### Experiment 2
#################################################################
###################################################################

### Clear workspace to remove first experiment
rm(list=ls())

### Set iterations

##### Libraries
library(Hmisc)
library(gtools)
library(gplots)
library(ggplot2)

#### Globals
nIter <- 5000
random.seed <- 1628 # Generated via "Random.org" [1-5000]

## Load in results file
data <- read.csv("experiment2_public.csv", stringsAsFactors=F)

## Non-response is coded as -99 for seen questions in raw data

#################################
## Data pre-processing
#################################

## Make respondent ID a character
data$ResponseId <- as.character(data$ResponseId)

## Merge responses
## Damages responses
data$damages <- data$damages1 
data$damages[data$damages1 == ""] <- data$damages2[data$damages1 == ""]
data$damages[data$damages == ""] <- NA

data$damagesClaim <- ifelse(grepl("Claimant", data$damages), 1, 0)
data$damagesResp <- ifelse(grepl("Respondent", data$damages), 1, 0)
data$damagesClaim[is.na(data$damages)] <- NA
data$damagesResp[is.na(data$damages)] <- NA

## Expenses Responses
data$expenses <- data$expenses1
data$expenses[data$expenses1 == ""] <- data$expenses2[data$expenses1 == ""]
data$expenses[data$expenses == ""] <- NA

data$expensesSplit <- ifelse(data$expenses == "Each party bears its own expenses.", 1, 0)
data$expensesSome <- ifelse(data$expenses == "Respondent reimburses some of the Claimant's expenses.", 1, 0)
data$expensesAll <- ifelse(data$expenses == "Respondent reimburses all of the Claimant's expenses.", 1, 0)
data$expensesWrong <- ifelse(grepl("Claimant reimburses", data$expenses), 1, 0)
data$expensesAny <- as.integer(data$expensesSome==1|data$expensesAll==1)

## Create indicators for appointer
data$AppointedResp <- ifelse(data$Appointer == " by the Respondent.", 1, 0)
data$AppointedClaim <- ifelse(data$Appointer == " by the Claimant.", 1, 0)
data$AppointedParties <- ifelse(data$Appointer == " by the Parties.", 1, 0)
data$AppointedBlind <- ifelse(data$Appointer == ".", 1, 0) ### Dot denotes blind appointment

data$Appointer <- relevel(factor(data$Appointer), " by the Claimant.")
data$Appointer <- relevel(factor(data$Appointer), ".") ### Dot denotes blind appointment

### Did they finish the survey?
data$completed <- ifelse(data$Finished == TRUE,1, 0)

#### Drop respondents that did not finish or said "no" to the consent form
data.finished = subset(data, Finished == TRUE&consent=="Yes")

### Subset out respondents that did not answer the outcome question.
data.finished.damages <- subset(data.finished, !is.na(data.finished$damagesClaim))
data.finished.expenses <- subset(data.finished, !is.na(data.finished$expenses))

#### How many observations in each category?
table(data.finished.expenses$Appointer)
table(data.finished.expenses$RespAward)

table(data.finished.damages$Appointer)
table(data.finished.damages$RespAward)

### Check for whether randomization holds, p > .10 in all cases.
chisq.test(table(data.finished.expenses$Appointer), p=c(1/6, 1/3, 1/6, 1/3)) ## Appointer
chisq.test(table(data.finished.expenses$RespAward), p=c(1/3, 1/3, 1/3)) ## Respondent's request

chisq.test(table(data.finished.damages$Appointer), p=c(1/6, 1/3, 1/6, 1/3)) ## Appointer
chisq.test(table(data.finished.damages$RespAward), p=c(1/3, 1/3, 1/3)) ## Respondent's request

###############################
### Bootstrap and make plots
###############################

########################################
#### Plot of marginal effects for Claimant v. Respondent case
########################################

#### Treatment Effects

### Placeholder frame to store effect estimates
results_frame_exp2 <- NULL

######### Bootstrap

## Set random seed
set.seed(random.seed)

## Matrix to store marginal effect
appointerRespondent_Damages <- matrix(nrow=nIter, ncol=1)

data.finished.damages$Appointer <- relevel(data.finished.damages$Appointer, " by the Respondent.")

# Bootstrap routine
for (iter in 1:nIter){

  # Resample rows of data
  if (iter != 1){ # Keep the data the same for the first estimate
    indices <- sample(1:nrow(data.finished.damages), nrow(data.finished.damages), replace=T) 
    boot.data <- data.finished.damages[indices,]
  }else{
    boot.data <- data.finished.damages
  }
  
  # Difference-in-means
  boot_damages <- lm(damagesClaim ~ Appointer, data=subset(boot.data))
  
  appointerRespondent_Damages[iter,] <- c(coefficients(boot_damages)["Appointer by the Claimant."])
}

## Point estimate
appointerRespondent_Damages_point <- appointerRespondent_Damages[1,]
## Standard Error
appointerRespondent_Damages_se <- apply(appointerRespondent_Damages, 2, function(x) sd(x))
## Confidence intervals
appointerRespondent_Damages_ci <- apply(appointerRespondent_Damages, 2, function(x) quantile(x, c(.025, .975)))
appointerRespondent_Damages_ci90 <- apply(appointerRespondent_Damages, 2, function(x) quantile(x, c(.05, .95)))
## P-values (asymptotic)
appointerRespondent_Damages_pval <- sapply(appointerRespondent_Damages_point/appointerRespondent_Damages_se, function(x) 2*pnorm(-abs(x)))


### Store in results frame
results_frame_exp2 <- rbind(results_frame_exp2,
      data.frame(treatment=c("Appointed by Claimant v. Respondent"), outcome=c("Damages (Pro-Claimant)"), estimate=round(appointerRespondent_Damages_point, 3)[1:1], 
                 standard_error = round(appointerRespondent_Damages_se, 3)[1:1], pval = round(appointerRespondent_Damages_pval, 3)[1:1]))

### Make a plot
damages_effect <- data.frame(choice = c( "Claimant's Damages\n(Baseline: Respondent's)"),
                             point= appointerRespondent_Damages_point,
                             lower95 = apply(appointerRespondent_Damages, 2, function(x) quantile(x, c(.025, .975)))[1,1], upper95=
                               apply(appointerRespondent_Damages, 2, function(x) quantile(x, c(.025, .975)))[2,1],
                             lower90 = apply(appointerRespondent_Damages, 2, function(x) quantile(x, c(.05, .95)))[1,1], upper90=
                               apply(appointerRespondent_Damages, 2, function(x) quantile(x, c(.05, .95)))[2,1])

png("appointed_respondent_damages_exp2.png", width=900, height=200)
p = ggplot(damages_effect, aes(y=point, x=choice))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower90,ymax=upper90), size= 4)
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Change in the probability respondents choose the\nclaimant's damages proposal", limits=c(-.4, .4))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + ggtitle("Appointed by the Claimant versus the Respondent")
p = p + theme_bw1()
print(p)
dev.off()
ggsave("appointed_respondent_damages_exp2.pdf", width=12, height=2.2857)
################ Appointer (claimant v. respondent) - on expenses

## Set random seed
set.seed(random.seed)

## Matrix to store marginal effect
appointerRespondent_Expenses <- matrix(nrow=nIter, ncol=4)

data.finished.expenses$Appointer <- relevel(data.finished.expenses$Appointer, " by the Respondent.")

# Bootstrap routine
for (iter in 1:nIter){

  # Resample rows of data
  if (iter != 1){ # Keep the data the same for the first estimate
    indices <- sample(1:nrow(data.finished.expenses), nrow(data.finished.expenses), replace=T) 
    boot.data <- data.finished.expenses[indices,]
  }else{
    boot.data <- data.finished.expenses
  }
  
  # Difference-in-means
  
  boot_split <- lm(expensesSplit ~ Appointer, data=boot.data)
  boot_some <- lm(expensesSome ~ Appointer, data=boot.data)
  boot_all <- lm(expensesAll ~ Appointer, data=boot.data)
  boot_any <- lm(expensesAny ~ Appointer, data=boot.data)
  
  appointerRespondent_Expenses[iter,] <- c(coefficients(boot_all)["Appointer by the Claimant."], coefficients(boot_some)["Appointer by the Claimant."], coefficients(boot_split)["Appointer by the Claimant."],
                                           coefficients(boot_any)["Appointer by the Claimant."])
}

# Point estimate
appointerRespondent_Expenses_point <- appointerRespondent_Expenses[1,]
# Standard error
appointerRespondent_Expenses_se <- apply(appointerRespondent_Expenses, 2, function(x) sd(x))
# Confidence intervals
appointerRespondent_Expenses_ci <- apply(appointerRespondent_Expenses, 2, function(x) quantile(x, c(.025, .975)))
appointerRespondent_Expenses_ci90 <- apply(appointerRespondent_Expenses, 2, function(x) quantile(x, c(.05, .95)))
# asymptotic p-value
appointerRespondent_Expenses_pval <- sapply(appointerRespondent_Expenses_point/appointerRespondent_Expenses_se, function(x) 2*pnorm(-abs(x)))


### Store in results frame
results_frame_exp2 <- rbind(results_frame_exp2,
                            data.frame(treatment=c("Appointed by Claimant v. Respondent","Appointed by Claimant v. Respondent","Appointed by Claimant v. Respondent"), 
                                       outcome=c("All Costs", "Some Costs", "Split Costs"), estimate=round(appointerRespondent_Expenses_point, 3)[1:3], 
                                       standard_error = round(appointerRespondent_Expenses_se, 3)[1:3], pval = round(appointerRespondent_Expenses_pval, 3)[1:3]))



### Make a plot

expenses_effect <- data.frame(choice = c( "Respondent Reimburses All","Respondent Reimburses Some","Parties Pay Own Costs"),
                              point= appointerRespondent_Expenses[1,c(1:3)],
                              lower95 = apply(appointerRespondent_Expenses, 2, function(x) quantile(x, c(.025, .975)))[1,c(1:3)], upper95=
                                apply(appointerRespondent_Expenses, 2, function(x) quantile(x, c(.025, .975)))[2,c(1:3)],
                              lower90 = apply(appointerRespondent_Expenses, 2, function(x) quantile(x, c(.05, .95)))[1,c(1:3)], upper90=
                                apply(appointerRespondent_Expenses, 2, function(x) quantile(x, c(.05, .95)))[2,c(1:3)])

expenses_effect$choice <- factor(expenses_effect$choice, c("Respondent Reimburses All","Respondent Reimburses Some","Parties Pay Own Costs"))

png("appointed_respondent_expenses_exp2.png", width=900, height=350)
p = ggplot(expenses_effect, aes(y=point, x=choice))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower90,ymax=upper90), size= 4)
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Change in the probability respondents select the outcome", limits=c(-.4, .4))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + ggtitle("Appointed by the Claimant versus the Respondent")
p = p + theme_bw1()
print(p)
dev.off()
ggsave("appointed_respondent_expenses_exp2.pdf", width=12, height=4)

########################################
#### Plot of marginal effects for Parties v. Blind case
########################################

#### Treatment Effects

## Set random seed
set.seed(random.seed)

## Matrix to store marginal effect
appointerParties_Damages <- matrix(nrow=nIter, ncol=1)

data.finished.damages$Appointer <- relevel(data.finished.damages$Appointer, ".")

# Bootstrap routine
for (iter in 1:nIter){

  # Resample rows of data
  if (iter != 1){ # Keep the data the same for the first estimate
    indices <- sample(1:nrow(data.finished.damages), nrow(data.finished.damages), replace=T) 
    boot.data <- data.finished.damages[indices,]
  }else{
    boot.data <- data.finished.damages
  }
  
  # Difference-in-means
  boot_damages <- lm(damagesClaim ~ Appointer, data=subset(boot.data))
  
  appointerParties_Damages[iter,] <- c(coefficients(boot_damages)["Appointer by the Parties."])
}

appointerParties_Damages_point <- appointerParties_Damages[1,]
appointerParties_Damages_se <- apply(appointerParties_Damages, 2, function(x) sd(x))
appointerParties_Damages_ci <- apply(appointerParties_Damages, 2, function(x) quantile(x, c(.025, .975)))
appointerParties_Damages_ci90 <- apply(appointerParties_Damages, 2, function(x) quantile(x, c(.05, .95)))
appointerParties_Damages_pval <- sapply(appointerParties_Damages_point/appointerParties_Damages_se, function(x) 2*pnorm(-abs(x)))

### Store in results frame
results_frame_exp2 <- rbind(results_frame_exp2,
                            data.frame(treatment=c("Appointed by Parties v. Blind"), outcome=c("Damages (Pro-Claimant)"), estimate=round(appointerParties_Damages_point, 3)[1:1], 
                                       standard_error = round(appointerParties_Damages_se, 3)[1:1], pval = round(appointerParties_Damages_pval, 3)[1:1]))


## Make a plot
damages_effect_parties <- data.frame(choice = c( "Claimant's Damages\n(Baseline: Respondent's)"),
                             point= appointerParties_Damages_point,
                             lower95 = apply(appointerParties_Damages, 2, function(x) quantile(x, c(.025, .975)))[1,1], upper95=
                               apply(appointerParties_Damages, 2, function(x) quantile(x, c(.025, .975)))[2,1],
                             lower90 = apply(appointerParties_Damages, 2, function(x) quantile(x, c(.05, .95)))[1,1], upper90=
                               apply(appointerParties_Damages, 2, function(x) quantile(x, c(.05, .95)))[2,1])


png("appointed_parties_damages_exp2.png", width=900, height=200)
p = ggplot(damages_effect_parties, aes(y=point, x=choice))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower90,ymax=upper90), size= 4)
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Change in the probability respondents choose the\nclaimant's damages proposal", limits=c(-.4, .4))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + ggtitle("Appointed by the Parties versus Blind")
p = p + theme_bw1()
print(p)
dev.off()
ggsave("appointed_parties_damages_exp2.pdf", width=12, height=2.2857)

########## Expenses (Parties v. Blind)

## Set random seed
set.seed(random.seed)

## Matrix to store marginal effect
appointerParties_Expenses <- matrix(nrow=nIter, ncol=4)

data.finished.expenses$Appointer <- relevel(data.finished.expenses$Appointer, ".")

# Bootstrap routine
for (iter in 1:nIter){

  # Resample rows of data
  if (iter != 1){ # Keep the data the same for the first estimate
    indices <- sample(1:nrow(data.finished.expenses), nrow(data.finished.expenses), replace=T) 
    boot.data <- data.finished.expenses[indices,]
  }else{
    boot.data <- data.finished.expenses
  }
  
  # Difference-in-means
  
  boot_split <- lm(expensesSplit ~ Appointer, data=boot.data)
  boot_some <- lm(expensesSome ~ Appointer, data=boot.data)
  boot_all <- lm(expensesAll ~ Appointer, data=boot.data)
  boot_any <- lm(expensesAny ~ Appointer, data=boot.data)
  
  appointerParties_Expenses[iter,] <- c(coefficients(boot_all)["Appointer by the Parties."], coefficients(boot_some)["Appointer by the Parties."], coefficients(boot_split)["Appointer by the Parties."],
                                           coefficients(boot_any)["Appointer by the Parties."])
}

appointerParties_Expenses_point <- appointerParties_Expenses[1,]
appointerParties_Expenses_se <- apply(appointerParties_Expenses, 2, function(x) sd(x))
appointerParties_Expenses_ci <- apply(appointerParties_Expenses, 2, function(x) quantile(x, c(.025, .975)))
appointerParties_Expenses_ci90 <- apply(appointerParties_Expenses, 2, function(x) quantile(x, c(.05, .95)))
appointerParties_Expenses_pval <- sapply(appointerParties_Expenses_point/appointerParties_Expenses_se, function(x) 2*pnorm(-abs(x)))

### Store in results frame
results_frame_exp2 <- rbind(results_frame_exp2,
                            data.frame(treatment=c("Appointed by Parties v. Blind","Appointed by Parties v. Blind","Appointed by Parties v. Blind"), 
                                       outcome=c("All Costs", "Some Costs", "Split Costs"), estimate=round(appointerParties_Expenses_point, 3)[1:3], 
                                       standard_error = round(appointerParties_Expenses_se, 3)[1:3], pval = round(appointerParties_Expenses_pval, 3)[1:3]))


### Make a plot
expenses_effect_parties <- data.frame(choice = c( "Respondent Reimburses All","Respondent Reimburses Some","Parties Pay Own Costs"),
                              point= appointerParties_Expenses[1,c(1:3)],
                              lower95 = apply(appointerParties_Expenses, 2, function(x) quantile(x, c(.025, .975)))[1,c(1:3)], upper95=
                                apply(appointerParties_Expenses, 2, function(x) quantile(x, c(.025, .975)))[2,c(1:3)],
                              lower90 = apply(appointerParties_Expenses, 2, function(x) quantile(x, c(.05, .95)))[1,c(1:3)], upper90=
                                apply(appointerParties_Expenses, 2, function(x) quantile(x, c(.05, .95)))[2,c(1:3)])

expenses_effect_parties$choice <- factor(expenses_effect_parties$choice, c("Respondent Reimburses All","Respondent Reimburses Some","Parties Pay Own Costs"))

png("appointed_parties_expenses_exp2.png", width=900, height=350)
p = ggplot(expenses_effect_parties, aes(y=point, x=choice))
p = p + coord_flip()  
p = p + geom_hline(yintercept = 0,size=1,linetype="dashed") 
p = p + geom_pointrange(aes(ymin=lower90,ymax=upper90), size= 4)
p = p + geom_pointrange(aes(ymin=lower95,ymax=upper95), size= 2)
p = p + scale_y_continuous(name="Change in the probability respondents select the outcome", limits=c(-.5, .5))
p = p + scale_x_discrete(name="")

theme_bw1 <- function(base_size = 20, base_family = "") {
  
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size*.9, colour = "black"),
      axis.text.y =       element_text(size = base_size, colour = "black"), # changes position of X axis text
      
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      plot.title =             element_text(face = "bold"),
      legend.position = "none"
    )
}
p = p + ggtitle("Appointed by the Parties versus Blind")
p = p + theme_bw1()
print(p)
dev.off()
ggsave("appointed_parties_expenses_exp2.pdf", width=12, height=4)

#### Finally, print point estimates, SEs and p-values for effect estimates in Experiment 2

print(results_frame_exp2)

### END EXPERIMENT 2